Minimizing chi-square by taking into account the error bars

Needs["ErrorBarPlots`"]

xi={0.8,2.0,3.0,3.8,5.1,6.1};
y={4.3,5.6,12.3,13.0,18.7,19.1};
σ={0.5,2.1,1.3,2.3,2.8,0.6};

data=Table[{xi[[i]],y[[i]]},{i,1,Last[xi]}];

ListPlot[data,AxesOrigin->{0,0},PlotStyle->PointSize[0.015]]

data2=Table[{data[[i]],ErrorBar[σ[[i]]]},{i,1,Length[σ]}]

p1=ErrorListPlot[data2,PlotStyle->PointSize[0.015],PlotRange->{0,25}]

param=FindFit[data,m x+b,{m,b},x]

p2=Plot[m x+b/.param,{x,0,7},PlotLabel->" m=3.113,    b=1.374"]

χsquare=(∑_(i=1)^(Length[xi]))(((((y[[i]]-(m xi[[i]]+b)/.param)/(σ[[i]]))))^2)

Here is the best fit when error bars are not included in the fitting procedure.

Show[{p1,p2},PlotLabel->"(χ^2)=7.87,   m=3.113,   b=1.374"]

Let’s formulat a (χ^2) function that takes into account both variables (m and b)

χ2[m_,b_]:=(∑_(i=1)^(Length[xi]))(((((y[[i]]-(m xi[[i]]+b))/(σ[[i]]))))^2)

Here is the (χ^2)(b) where m has the value from FindFit from up above.
This does not represent the “best fit” to “b”  but we’ll investigate it nonetheless.

χ2b[b_]:=(∑_(i=1)^(Length[xi]))(((((y[[i]]-((m/.param ) xi[[i]]+b))/(σ[[i]]))))^2)

g1=Plot[χ2b[b],{b,0.3,2.3},PlotRange->{6,11},PlotLabel->"(χ^2)  vs.  b (intercept)"];

Here is the (χ^2)(m) where b has the value from FindFit from up above.
This does not represent the “best fit”  to “m” but we’ll investigate it nonetheless.

χ2m[m_]:=(∑_(i=1)^(Length[xi]))(((((y[[i]]-(m xi[[i]]+(b/.param)))/(σ[[i]]))))^2)

g2=Plot[χ2m[m],{m,2.7,3.2},PlotRange->{4.0,8.0},PlotLabel->"(χ^2)  vs.  m (slope)"];

Note:  The (χ^(2 ))'s do not share the same minimum.  That’s because the parameters “a” and “b” were not calculated with the error bar information.

GraphicsRow[{g1,g2}]

Let’s minimize the (χ^2) by varying the slope (m) and intercept (b) simultaneously.
This represents the best approach for finding the optimum parameters for the slope (m) and the intercept (b).  FindMinimum is a Mathematica function that finds the minimum of a multivariate function.  In the example below, FindMinimum uses the “m” and “b” values from FindFit (up above) as a starting point for the fit.

param2=FindMinimum[χ2[m,b],{{m,m/.param},{b,b/.param}}]

param2[[2]]

p3=Plot[(m x+b)/.param2[[2]],{x,0,7},PlotLabel->"(χ^2)=3.105,  m=2.932,  b=2.124"]

g4=Show[{p1,p3},PlotLabel->"(χ^2)=3.589,  m=2.817,  b=2.119      (using error-bars)"];

g3=Show[{p1,p2},PlotLabel->"(χ^2)=7.87,  m=3.113,  b=1.374     (not using error-bars)"];

Here is a comparison between the two fits:

On the left:      The fit was made without using error-bars while using the FindFit function
		 The (χ^2) was calculated using the error-bars.
On the right:	 The fit was made using error-bars while using the FindMinimum function
		 The (χ^2) was calculated using the error-bars.

GraphicsRow[{g3,g4}]

On the left:	Since the error-bars were not used in this fit, the (χ^2)= 7.87 does
		not have much meaning.  There is no useful interpretation for this (χ^2).
On the right:	Error-bars were used in this fit, so, the (χ^2) and reduced-(χ^2) are both 
		meaningful statistical quantites.  The reduced-(χ^2) in this case is the (χ^2) 
		divided by the degrees of freedom (6-2=4).  As a result, the 
		reduced-(χ^2)=0.897 is a meaningful statistical quantity.
		
Note:		The parameters “m” and “b” are noticeably different between the two
		fits.  This is another feature that occurs when using the error-bars 
		correctly.

Let’s perform the error analysis using the correct values for “m” and “b.”  Using the values from the FindMinimum fit we find:

χ2bFM[b_]:=(∑_(i=1)^(Length[xi]))(((((y[[i]]-((m/.param2[[2]] ) xi[[i]]+b))/(σ[[i]]))))^2)
χ2mFM[m_]:=(∑_(i=1)^(Length[xi]))(((((y[[i]]-(m xi[[i]]+(b/.param2[[2]])))/(σ[[i]]))))^2)
g5=Plot[χ2bFM[b],{b,1,3},PlotRange->{3,6},PlotLabel->"(χ^2)  vs.  b (intercept)"];
g6=Plot[χ2mFM[m],{m,2.5,3.1},PlotRange->{3,6},PlotLabel->"(χ^2)  vs.  m (slope)"];
GraphicsRow[{g5,g6}]

Note:	Both (χ^2) curves have the same minimum value, whereas before, they did not.

The errors in the slope and the intercept (δm and δb) can be found by one of the following techniques:

1.    	Use the determinants program and calculate δm and δb
	directly, since this is a linear fit.
	
2.	Use the (χ^2)+1  technique, and you should get the same result.
	Since technique #1 is only useful for linear fits, you must use
	a technique similar to the (χ^2)+1  technique if your fitting
	function is non-linear.

Let’s use technique #1 to find the errors in the slope in the intercept (δm and δb).

num=Length[xi];
delta =(∑_(i=1)^num)(1/((σ[[i]])^2)) (∑_(i=1)^num)(((xi[[i]])^2)/((σ[[i]])^2)) - ((((∑_(i=1)^num)((xi[[i]])/((σ[[i]])^2))))^2);
a1 = (1/delta)((∑_(i=1)^num)(((xi[[i]])^2)/((σ[[i]])^2)) (∑_(i=1)^num)((y[[i]])/((σ[[i]])^2)) - (∑_(i=1)^num)((xi[[i]])/((σ[[i]])^2)) (∑_(i=1)^num)((xi[[i]] y[[i]])/((σ[[i]])^2)));
b1 = (1/delta)((∑_(i=1)^num)(1/((σ[[i]])^2)) (∑_(i=1)^num)((xi[[i]] y[[i]])/((σ[[i]])^2)) - (∑_(i=1)^num)((xi[[i]])/((σ[[i]])^2)) (∑_(i=1)^num)((y[[i]])/((σ[[i]])^2)));
deltaa = (Sqrt[(1/delta)(∑_(i=1)^num)(((xi[[i]])^2)/((σ[[i]])^2))]);
deltab = (Sqrt[(1/delta)(∑_(i=1)^num)(1/((σ[[i]])^2))]);
Print["Intercept  a =  ",a1, "  ± ",deltaa]
Print["Slope      b = ",b1, " ± ",deltab]

Note:  the slopes and intercepts are identically the same when comparing the results between FindMinimum and the Determinant techniques.  The Determinant technique also gives the correct errors for the slope and the intercept (δm and δb).

If this were a non-linear fit, you need to use the FindMinimum approach and calculate the errors for the slope and intercept (δm and δb) using the 2nd technique above.

The lessons to be learned:

1.	Always include your error-bars when determining the values
	of your fitting parameters (e.g., m and b).
	
2.	Use either technique #1 or #2 described above to calculate
	the uncertainties of your fitting parameters (e.g., δm and δb).
	Technique #1 is easier to use if you’re performing a linear fit.